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Abstract 

We present an analytical study of the spatial decay 
rate 7 of the one-particle density matrix p(r, r') ^ 



exp(— 7|r - 



for systems described by single parti- 



cle orbitals in periodic potentials in arbitrary dimen- 
sions. This decay reflects electronic locality in con- 
densed matter systems and is also crucial for 0{N) 
density functional methods. We find that 7 behaves 
contrary to the conventional wisdom that generically 
7 cx -s/A in insulators and 7 cx \/T in metals, where 
A is the direct band gap and T the temperature. 
Rather, in semiconductors 7 cx A, and in metals at 
low temperature 7 oc T. 



Density functional theory (DFT)||i[ describes 
many-body systems via a single-particle formalism 
and is the basis for modern, large-scale calculations in 
solid-state systems [Q. The one-particle density ma- 
trix p = IV'n)/n(V'n|, which describes the state of 
a single-particle quantum system, is the key quantity 
needed for the computation of physical observables: 
total system energies, atomic forces, and phonons can 
all be computed directly from p. 

Remarkably, despite the de-localized nature of the 
single particle states |'0n), which may extend across 
an entire solid, the physics of the electronic states in 
a given region of a material is affected only by the 
local environment. Reflecting this, the force on an 
atom depends mostly on the positions of its nearest 
neighbors. This electronic localization is manifest in 
the "nearsightedness"!^ of p: p{r,f') = {f\p\f') ^ 
exp(— 7|r — f'\) where 7 > 0. This exponential decay 
has been verified numerically |j and analytically ||, 

The locality of p not only is important for under- 
standing the nearsightedness of effects arising from 
electronic structure but also has direct practical im- 



pact on DFT calculations. Recently, methods have 
been proposed 0, @ |§, |ll g H, H ||, g 
that use p directly and exploit its locality. Compu- 
tationally, these methods scale as 0{N), where N is 
the number of atoms in the simulation cell. However, 
their prefactors depend strongly on 7: some scale as 
iV/76|l|, 1^, H and others as N/-/^^. Knowing 
how 7 depends on the system under study is thus 
critical for carrying out such calculations. For a re- 
view of 0{N) methods, see fl^ . 

Generally, solid-state systems have an underlying 
periodic structure. The introduction of localized 
defects]^ or surfaces |9| does not change the spatial 
range of p from that of the underlying periodic lat- 
tice. Thus, understanding the locality of p even for 
perfectly periodic systems is of direct relevance for 
realistic material studies. To date, the generic be- 
havior of 7 is poorly understood. For insulators in 
one dimension, Kohn has shown that 7 oc \J —E^ in 
the tight-binding limit where is an atomic ion- 
ization energy p. Motivated by this, it has been 
assumed ^, 24 and argued that 7 cx \fK in 
multiple dimensions and more general conditions, 
where A is the band gap. For metals, it has been 
assumed [H and argued jis) that 7 cx -s/T, where T 
is the electronic temperature. However, the results 
in|p^, which to date represent the only effort to de- 
termine 7 generically, are based on the assumption 
that the inverse of the overlap matrix of a set of Gaus- 
sian orbitals decays in a Gaussian manner. On the 
contrary, the inverses of such overlap matrices decay 
only exponentially, and thus the behavior of 7 war- 
rants further study. 

Here, we show that the behavior of 7 is more com- 
plex than previously assumed. For insulators, 7 is 
determined by the analytical behavior of the filled 
bands, which is determined by the strength of the 
periodic potential which, in turn, is most strongly re- 
flected in the size of the direct band gaps. Indirect 
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gaps, being more accidentally related to the strength 
of the potential, have a more haphazard relation to 
7, which we do not consider here. 

For insulating systems, we find that 7 has the fol- 
lowing asymptotic behavior as a function of the direct 
gap A, lattice constant a, and electron mass m. 



p{r, r') are related via 



7 



111 



for a^A 
for a^A 







(weak-binding) 
(tight-binding) 



The indeterminacy in the tight-binding limit re- 
sults from the electronic states becoming atomic or- 
bitals, and thus, 7 depends on the details of the un- 
derlying atomic potential. Some systems (see below) 
exhibit 7 oc a/A, but this is not universal, as previ- 
ously assumed. 

One can, however, make a definitive statement in 
the weak-binding limit, which is of direct impor- 
tance for small-gap systems such as those with weak 
pseudopotentials or gaps due to Jahn- Teller distor- 
tions. For example, semiconductors such as Si and 
GaAs have gaps that are significantly smaller than 
their band widths, and we expect them to fall into 
the weak-binding case. In Si and GaAs, we find 
a^Am/?i^ 4 and ^ 2.5, respectively. Inspecting 
Figure ^ we see that for such values the behavior of 
7 is well in the weakly-bound limit. 

For metals with a fixed number of electrons, we 
find 



^Tln 



(^) 



for T 
for T 







(quantum) 
(classical) 



in terms of the temperature T, the typical gradient 
of the band energy eg on the Fermi surface, and the 
Fermi energy sp. 

The low-temperature result is of direct practical in- 
terest for calculations in metals. (This result was also 
found in a contemporaneously submitted publication 
[ 141.) We find behavior resembling that proposed in 
[4 |l5|, i.e. 7 (X Vr, only at extremely high temper- 
atures. 

We now present analytical arguments that substan- 
tiate the above results and shed light on the physical 
mechanisms leading to and differentiating among the 
different limits. We consider periodic systems with 
lattice vectors of characteristic length a. We choose 
units such that /m — 1 and kg = 1 in order to 
avoid cumbersome mathematical expressions; the re- 
sults presented above are easily recovered by insert- 
ing h^/m and fc^, as appropriate, in each step of the 
analysis below. 

The Bloch wave-functions V'„fei Wannier functions 
Wn, Fermi-Dirac fillings /^^r, and density matrix 



dk ( 



-ik-R 



Fn{R) 

Pn{r,f') 



^ J2 Wn (r, R)FniR~R') (r^ , R') 



R R' 



P{r,r') = ^p„(f,f'). 



(1) 



The integrals are over the first Brillouin zone with 
volume fls- R ranges over the lattice vectors. 
Pn{r,r') is the density matrix of the nth band, and 
p is a simple sum over all p„. Thus, we need only 
study the behavior of p„ for a given n. We analyze 
the behavior of 7 for the two possible cases of practi- 
cal interest, insulators at low temperature and metals 
at non-zero temperature. 

Insulators (T = 0) — When the chemical potential 
fi falls in the energy gap, all fillings are 1 or 0. A 
filled band with /^r = 1 has F{R) = q so that p is 
simply 

p(f,r") ^^W{f,R)W*if",R). 

R 

Wannier functions are exponentially localized |[ 
|, I, m |l|, |l| and satisfy VF(f, i? ) = W{f- R, 0). 
Thus, only a finite set of R contribute significantly to 
the above sum, and the decay rates of p and W are 
the same. Therefore, we need only determine 7 for 
the Wannier functions. 

As a concrete example and an initial orientation, 
we solve exactly for 7 for the lowest band of a model 
one dimensional system for all binding strengths. We 
choose the periodic potential to be that of an ar- 
ray of attractive delta-functions of strength V > 0, 
U{x) = —V^j^S{x — no). Following]^, we define 
p{e) = cos(a\/te) - 1/ sin(oV2e)/ V2e. The band 
structure is found by solving cos(fca) = pi^k) for 
real k. Focusing on the lowest band, we denote e as 
the value of e where p{e) achieves its first minimum. 
Kohn has shown that 7 = cosh^^ Im(£)I- We solve 
the above transcendental system numerically for dif- 
ferent values of V and plot as a function of a^A 
in Figure 0[a]. The behavior at small A is clearly 
linear, showing that 7 ~ aA for a weak potential. 
The leading asymptotic behavior is 7 '--^ \/A for a 
strong potential (i.e. large A). In this case, the gen- 
eral notion that 7 cx VA is clearly incorrect for weak 
potentials, and it is natural to ask whether this re- 
sult is peculiar to our simple model or whether it is 
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Figure 1: [a] 07 versus a^A for a periodic array of 
attractive delta potentials (solid curve). The dashed 
line is 07 = a^A/(27r) . The dotted curve is 07 = 
aVSA, the leading asymptotic behavior of 7 for large 
a^A. [b] a7 versus a^A for a cubic lattice of Gaussian 
potentials: 7 in the [100] (circles), [110] (stars) and 
[111] (pluses) directions. 

more universal. As we now argue, for a weak poten- 
tial, 7 aA is quite general, whereas in the case of a 
strong potential, the behavior of 7 is not unique and 
depends on the details of the atomic system underly- 
ing the periodic lattice. The crossover from weak to 
strong potential behavior should occur when A is of 
order of the band width. In the figure, this occurs for 
A ~ 5(7r/a)^. We now analyze each case separately. 

Weak-binding insulators in general — We wish to 
find 7 in the limit of a weak periodic potential U{r). 
Eqs. (|^) show that the Wannicr function is the 
Fourier transform of '(/'£• Thus the range 6k in k- 
space where ipj: has its strongest variations deter- 
mines the spatial range of W. From basic Fourier 
analysis, ^ ^ 6k. 

A simple heuristic argument shows that 7 ~ a A. 
Starting with a free electron description, a weak po- 
tential U{f) causes the opening of a gap A at the 
edges of the Brillouin zone. The extent 5k of the 
region about the zone-edges where deviates most 
appreciably from its free electron value is given by 
5k^ /2m* ^ A, where m* is the effective mass at the 
zone edge. Standard treatments [^ show that for 
weak potentials m* ~ a^A. Combining these results, 
we see that 6k ^ aA, whence 7 ~ aA. 

This heuristic argument gives the desired result, 
but there are hidden assumptions. The argument is 
based solely on the behavior of the band structure e^:. 



whereas W is determined by the wavefunctions ^g. 
One must be sure that £g and -0^ vary over the same 
range 6k. Thus, We present a more precise argument 
in terms of -0^: alone. 

Letting = e''^'^ u^{r) , Eqs. (|^) show that 

W is also a Fourier transform of u-^. Away from the 
edges of the Brillouin zone, is given perturbatively 

by 

^^[k^-\k + G\^] /2 

where {r\G) — e*'^ '" and G is a reciprocal lattice 
vector: uj: is smooth and analytic in fc, and uj: ^ 1. 

However, close to the zone edges, « |fc-l-G| and 
deviates appreciably from unity in a region satisfying 
[fc2 - \k + G|2]/2 < V, where V is the typical size 
of the matrix elements of U. Since |G| ~ 1/a, this 
region has a width 6k ~ aV ^ a A and hence 7 ~ aA. 

As a concrete example, we study a cubic lattice of 
attractive Gaussian potentials with rms width a/n. 
We vary the depth of the potential, and for each 
depth, we compute A and p by sampling the Bril- 
louin zone on a cubic grid of size 40'^ and expanding 
ipj: in plane waves with |G| < 12-^. Diagonalizing 
the resulting Hamiltonian gives the ground-state ipj: 
from which we compute the density matrix. Sam- 
pling p(0,f') in the [100], [110] and [111] directions 
gives exponentially decaying envelopes upon which 
we perform linear fits on log plots to extract 7. Fig- 
ure Qb] shows our results, from which the behavior 
7 oc A is evident. 

Tight-binding insulators in general — The poten- 
tial U is the periodic sum of an atomic potential Vat, 
U{r) = J^R^atir- R). For sufficiently strong Vat, 
system properties are determined by the atomic po- 
tential. The Wannier functions become atomic or- 
bitals localized about the minima of Vat- Now, 7 
depends on the details of Vat and no single universal 
scaling can be found. To demonstrate the complexity 
and richness of this limit, we discuss briefly different 
examples of atomic potentials that lead to differing 
forms for 7. Note that in this atomic limit the lattice 
constant a is irrelevant in determining 7. 

For the Coulomb potential, Vat(r) — —Ze^jr. 
In the limit Ze^ — > 00, we have hydrogenic states 
centered on the lattice sites with energies i?„ = 
—Z'^e^ I2v? and Bohr radii = n? /Ze^. The gap 
A is an energy difference between atomic states, and 
so A ^ Z^e'^. Also, 7 ~ clq^ , and so we conclude 
7 ~ a/a. More generally, for any atomic potential 
with only a single dimensionful parameter (e.g. Ze^ 
above), dimensional analysis gives 7 ~ \/A. 
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Figure 2: 07 versus aTrTiJ,~2 for a BCC lattice of 
tight-binding s-orbitals: 7 in the [100] (circles), [110] 
(stars) and [111] (pluses) directions, /i is measured 
from the bottom of the band (see text). 



However, a similar analysis applied to a Gaus- 
sian potential, Vat{r ) = —Ve~^ 1"^" , gives 7 ^ Act, 
whereas, for a spherical well, Vat(r ) = —V6{a — r), 
we find that 7 has no dependence on A. Thus in 
the tight-binding limit, it is difficult to make generic 
statements regarding 7. 

Metals (T > Q) — In metals, the fillings ex- 
hibit rapid variations in k across the Fermi surface 
and F[R ) in (Q) becomes long-ranged. The Wannier 
functions, being independent of the fillings, remain 
exponentially localized, as discussed above. These 
facts combined with the structure of the sum in (0) 
imply that 7 in this case is determined by F[R). 
For an initial orientation, consider a band with a 
free electron-like form Sk = whose spherical 

Fermi surface is contained inside the first Brillouin 
zone. For such a band, F{R) is given by (|^) with 
/fc = 1/ (1 + exp [(fcV2 - ^J.) /T] ) . Below, we will use 
the fact that the density matrix of a true free electron 
gas is proportional to this F: p(f, r') oc F(f — r '). 

Because the Fermi surface is contained within the 
first zone, we may extend the k integral for F to in- 
finity. Changing to spherical coordinates, integrating 
by parts and using trigonometric identities yields 



f{r: 



1 



RdR 



dk cos (kR) 



00 cosh 



2T 



(2) 

When closing the integral in the upper complex k- 
plane, the relevant poles of the integrand are at fc = 
ki = ±^j2^i±2inT{2l+ 1) for integers I > 0. The 

residues of these poles contain the factor e''^'^ which 



gives rise to oscillations due to the real part of ki and 
exponential decay due to its imaginary part. 

When T — !■ 0, /i equals the Fermi energy fi — ep — 
kp/2 where kp = {Sn^ny^^ and n is the electron 
density. Thus we have ki « mT{2l -I- l)/kp ± kp- As 
R 00, the / = contribution dominates, so that 
7 = TrT/kp. 

As T ^ 00, the ideal gas result fi w Tln{nX^) 
holds where At = ^/2tt/T is the thermal de Broglie 
wavelength. In this limit, ki w so that 7 = 

/m^2rin(nA^) - y/TlTi{T/ep). 

Note that if we approximate the integrand of 
Eq. ^ by e^^/'^~'^ which corresponds to using 
Maxwell-Boltzmann fillings, F will have a Gaussian 
form F{R) oc e~^^ However, one can show that 
this approximation is only valid for small R. For 
large R, an exponential tail e^^^ remains where 7 is 
as described above. 

We now present separate arguments showing that 
these asymptotic forms for 7 are correct for metals in 
general. 

Metals asT — We first consider first T = 0. 
Bands below the Fermi level then have = 1 and 
F{R) = S^Q, and, as for the insulating case, their 
density matrices decay exponentially. However, for 
bands that cross the Fermi level, jumps discontin- 
uously from unity to zero wherever = fj,. As is well 
known, the Fourier transform of a discontinuous func- 
tion has algebraic falloff, and thus F{R)\t=o oc |i? I"** 
where 77 > 0. Such bands therefore dominate the de- 
cay of p as T — > 0. 

At finite T, the fillings are f^={l + e^)~^ where 
y = (e^ — ji) /T. As T ^ 0, the fillings now go from 
unity to zero in a narrow region about the Fermi sur- 
face defined by vectors kp satisfying \e^ — ~ T. To 
determine the width of this region, we approximate 
about the Fermi vector kp via y ~ Ve • (fc — kp)/T. 
The width of the transition region and 7 thus are 
given by 7 ~ (5fc ^ r/|Ve|. This argument holds for 
any Fermi surface no matter how complex (metallic 
or semi-metallic) at sufficiently low T such that 5k is 
smaller than the typical scale of the features of the 
Fermi surface. 

As a further verification for the general case, we 
study a body-centered cubic lattice of tight-binding 
s-orbitals with lattice constant a = 4.32A, for which 
the Fermi surface is non-spherical. We choose the 
tight-binding matrix element so that the band struc- 
ture has the free electron effective mass at fc = 0. 
Choosing to be 2/5 of the way from the band min- 
imum to the band maximum, we calculate F{R ) for 
various values of T . Sampling the Brillouin zone on 
a 200"^ grid, plotting F{R ) and finding exponentially 
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decaying envelopes, we perform linear fits on a log 
plot and extract 7. Figure || shows that indeed 7 (x T 
for T^O. 

Metals as T —> 00 — Here the electron kinetic en- 
ergy is much larger than the periodic potential so 
that we may approximate eg = the system is a 

classical ideal gas and is not of interest for solid-state 
calculations. Our previously derived result for a free 
electron gas yields 7 ^ ln(T/e_F)- Only in this 
limit do we find a result resembling that of . 
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